Make your own free website on Tripod.com


 algorithms computer graphics maths technical documents 




Algorithms





What follows is a chaotic collection of some maths and explanations of mainly (say 98 %) computer graphics algorithms....
Most stuff are my personal saved items from the comp.graphics.algorithms newsgroup. See what you can do with it
.






.... 3d rotations ... !!Gimbal lock?!! Don't let it get you down! In the following text, we will take the example of transforming the points data of a spaceship. We've got our spaceship aligned down our z. First of all, we want to roll the spaceship... If we do it now, it saves about 5000% calculation. Thats serious cpu. So... we do our standard rotate about the Z... I'm SURE you've got the calculation for that... But hell, we're all in the same boat, here it is: This is a matrix. The lines are to help me. +-------+-------+-------+ | cos t |-sin t | 0 | +-------+-------+-------+ | sin t | cos t | 0 | +-------+-------+-------+ | 0 | 0 | 1 | +-------+-------+-------+ Now, being a sane rational person, you'll look at that, see all the 0 and 1, and notice that you have four multiplies to do. Remember that: '4' four (4) This comes in later This method implies that our roll is an angle... not a bad assumption at all. Now, we have some vector somewhere which is the orientation of the plane.. the trickiest bit we'll encounter here, is actually calculating the axis... We're REALLY gonna need a LOT of bits in fix to store this, 'cos this is the most (in fact the only thing) prone to deterioration... We take this to be some axis that the user can control... The real hard work starts because we have to use rotate this by whatever, and it'll be quite complex to do interactively... However, suppose we do have this vector, as we almost certainly will have... We can take this vector, and call it 'd', because its the [d]irection the plane is pointing... Now, we want to map the z onto d... So, we rotate about our arbitrary axis... Uh oh! First of all, we must calculate the cross product of the two vectors... Uh oh! six multiplies! But... its only a preprocess step.. HOLD it right there! ALL of those six multiplies will be by either 0 or 1 Remember, Z is always {0,0,1} So... Lets quote the cross product: aXb= ay*bz - by*az bx*az - ax*bz ax*by - bx*ay Now, lets say that b is the z... So, we can insert the values to yield: aXz = ay*1 - 0*az 0*az - ax*1 ax*0 - 0*ay Well, this quickly simplifies to: aXz = ay -ax 0 This makes things a LOT easier.. Two variable assignments, and the whole cross product is done... Now, we've preprocessed this, we've still got that arbitrary axis rotation to do... But how many degrees to rotate it around?!??!?!? Well, we know that: A.B = |A||B|cos angle_between A and B are both unit vectors, so: A.B=cos angle between... But B=z So... az=cos angle between You can quickly verify this for yourself... So, to recap: We want to rotate about {ay,-ax,0} by angle arc-cosine(az) Now... Here's the matrix for a rotation about an arbitrary axis... I calculated this with a little help from a maths tutor, who found some info on quarternions in "the mathematical gazette" (don't ask ;-) and then simplified it into english.... I found out how to express an axis in quarternions, and then how to rotate a quarternion... The inns and outs were rather complex, and involved seven A4 sides of working... The things I will do for fun. Anyway, it ended up as a huge expansion, which concurs with Foley van Dam, Feiner, Hughes, et al... I shall quote theirs, as it uses some trigonometric identities, which make it more eloquent: We'll say our axis={x,y,z} for typing's sake I just hope my squared symbol shows up... Just a check: x = x*x that work? Also, t=theta Matrix: x+cos t(1-x) xy(1-cos t)-zsin t xz(1-cos t)+ysin t xy(1-cos t)+zsin t y+cos t(1-y) yz(1-cos t)-xsin t xz(1-cos t)-ysin t yz(1-cos t)+xsin t z+cos t(1-z) Hopefully, you can make sense of that... However, we can simplify it. Remember that the z component of our arbitrary axis was 0... That means that z in the above = 0... Well then, this yields: x+cos t(1-x) xy(1-cos t) ysin t xy(1-cos t) y+cos t(1-y) xsin t -ysin t xsin t cos t (I have merely overwritten the *z parts with ' ') Straightening this up yields: x+cos t(1-x) xy(1-cos t) ysin t xy(1-cos t) y+cos t(1-y) xsin t -ysin t xsin t cos t Now, this is very neat! This gives us nine multiplies per point... add that nine to the four from above... eleven multiplies to produce an aircraft at any orientation! Some people (me especially) hate stuff presented in matrix form, so i'll multiply out: for point {a,b,c} and axis {x,y,z} and angle t a'=a(x+cost(1-x)) + b(xy(1-cos t)) - c(ysin t) b'=a(xy(1-cost)) + b(y+cost(1-y)) + c(xsin t) c'=a(ysin t) - b(xsin t) + c(cos t) Et voila! However, you will note that we had two matrix form transformations, one for the roll, and another for the arbitrary. It +is+ possible to multiply these, but I haven't, just so that you can derive some mathematical enjoyment from this "experience" (also, i can't be bothered to multiply them :-) Source: I've provided some sample source. This is ONLY so that you can get a further idea of how to implement, it's relatively simple, and optimised towards simplifying the instructions used... It just offers a couple of functions... Don't expect this to work... It might, it might not... (it's not actually part of my implementation, as that was in ARM assembler - little use to most ppl here!) So, here it is: [ // +-------------------------------+ // | SpaceShip transformation code | // | Code to rotate a 3d object | // | with 'behaviour' similar to | // | a spaceship or aircraft | // +-------------------------------+ // This code specifies that the spaceship has been designed // as relative to the z axis... // i.e. that the point it rotates around is at the origin, // and its nose lies on the z axis. /* Y -Z | / /\|/____X \ \ / /\/ |_/ Z Sorry, but its hard (for me) to do a spaceship in ascii! */ // Usage of functions: // void Preprocess_CalculateMatrix(float x,float y,float z) // Use this to calculate the matrix needed to transform all the points // pass it the NORMALISED components of the vector describing // the direction of the nose from the origin. // (i.e. to point it straight towards you, send {0,0,1}) // void SetRoll(float angle) // Use this to set up the angle of roll (in radians, because thats what math.h wants) // Point Rotate_This_Point(Point) // use this to rotate a point. // I've just 'invented' a type here typedef float Point[3]; // this is just to help with parameter passing... // you don't need to use ANY of these variables, they're just for internal working float Rot_Matrix[3][3]; float ax,ay; float st,ct; float zs,zc; float x,y,z; // sorry these aren't very descriptive, but they get used a lot, // and it looks horrible otherwise void Preprocess_CalculateMatrix(float x, float y, float z) { // right then... // first off, sorry for using c++ comments // secondarily, the input values: // x,y,z are the *****NORMALIZED****** components of the vector // pointing in the direction which you want your object to now be aligned // IMPORTANT: Initially, the object is aligned down the Z axis ax=y; ay=-x; ct=z; st=sin(acs(ct)); // WARNING this line seems to imply a limitation that the this cannot face // backwards... there ought to be some form of test which would fix this, // but my trig seems to have gone to sleep. // anyway... not too hard to invert sign of all z components ;-) Rot_Matrix[0][0]=((ax*ax)+(ct*(1-(ax*ax)))); Rot_Matrix[1][0]=((ax*ay*(1-ct))); Rot_Matrix[2][0]=(ay*st); Rot_Matrix[0][1]=((ax*ay*(1-ct))); Rot_Matrix[1][1]=((ay*ay)+(ct*(1-(ay*ay)))); Rot_Matrix[2][1]=-(ax*st); Rot_Matrix[0][2]=-(ay*st); Rot_Matrix[1][2]=(ax*st); Rot_Matrix[2][2]=(ct); } void SetRoll(float angle) { zs=sin(angle); zc=cos(angle); } Point RotatePoint(Point p) { x=p[0]; y=p[1]; z=p[2]; // roll p[0]=(zc*x)+(zs*y); p[1]=(zc*y)-(zs*x); x=p[0]; y=p[1]; p[0]=(x*Rot_Matrix[0][0])+(y*Rot_Matrix[0][1])+(z*Rot_Matrix[0][2]); p[1]=(x*Rot_Matrix[1][0])+(y*Rot_Matrix[1][1])+(z*Rot_Matrix[1][2]); p[2]=(x*Rot_Matrix[2][0])+(y*Rot_Matrix[2][1])+(z*Rot_Matrix[2][2]); return p; } // and thats it... ] Again, this is just a (small) extension to my code, to give you a little idea of how this could be implemented... I offer no guarantee that this will work!








I have made a small proggy just to check out how I could do a blur. I used xor to get a nice looking image, and then calculated the average-color of the 4 surrounding colors. Is this the right way to do it ? It seems to go slow as hell... Can anyone give me any tip on how to increase the speed of the blur ? (in pascal or somewhat readable asm...) Anders Pedersen Code posted under: program _xorblur; uses crt,dos; TYPE Virtual = Array [1..64000] of byte; { The size of our Virtual Screen } VirtPtr = ^Virtual; { Pointer to the virtual screen } var x,y,i:integer; col:byte; Virscr,virscr2 : VirtPtr; { Our first Virtual screen } Vaddr,vaddr2 : word; { The segment of our virtual screen} procedure flip32(source,dest:Word); assembler; { This copies the entire screen at "source" to destination } asm push ds mov ax, [Dest] mov es, ax mov ax, [Source] mov ds, ax xor si, si xor di, di mov cx, 16000 db $F3, $66, $A5 pop ds end; Procedure Cls32 (Where:word;Col : Byte); assembler; { This clears the screen to the specified color } asm push es mov cx, 16000; mov es,[where] xor di,di mov al,[col] mov ah,al mov dx, ax db $66, $C1, $E0, $10 {shl eax, 16} mov ax, dx db $F3, $66, $AB {rep stosd} pop es End; Procedure SetUpVirtual; { This sets up the memory needed for the virtual screen } BEGIN GetMem (VirScr,64000); vaddr := seg (virscr^); {The segment for page 1} GetMem (VirScr2,64000); vaddr2 := seg (virscr2^); {The segment for page 2} END; {Setupvirtual} Procedure ShutDown; { This frees the memory used by the virtual screen } BEGIN FreeMem (VirScr,64000); FreeMem (VirScr2,64000); END; {Shutdown} procedure init;assembler; asm mov ax,13h int 10h end; procedure setpal(nr,r,g,b:byte); begin; port[$3c8]:=nr; port[$3c9]:=r; port[$3c9]:=g; port[$3c9]:=b; end; procedure putpixel(x,y:integer;col:byte); begin; mem[vaddr:x+(y*320)]:=col; end; procedure blur; var c,c1,c2,c3,c4:integer; color:single; begin; for y:=1 to 191 do for x:=1 to 320 do begin; c:=mem[$A000:x+(320*y)]; c1:=mem[$A000:(x-1)+(320*y)]; c2:=mem[$A000:(x+1)+(320*y)]; c3:=mem[$A000:x+(320*(y-1))]; c4:=mem[$A000:x+(320*(y+1))]; color:=round((c1+c2+c3+c4)/4); mem[vaddr:x+(320*y)]:=round(color); end; end; begin; init; setupvirtual; for i:=1 to 256 do setpal(i,0,i,i div 2); cls32(vaddr,0); for y:=1 to 191 do for x:=1 to 319 do putpixel(x,y,(x xor y)); flip32(vaddr,sega000); repeat; cls32(vaddr,0); blur; flip32(vaddr,sega000); until keypressed; shutdown; asm mov ax,3h;int 10h;end; end.










I my ray tracer I transform the direction of the starting ray like this: Vector vpn( Normalize( world->camera->target - world->camera->position ) ); Vector vup( 0,-1,0 ); vup = Normalize( vup- vpn*DotProd( vup, vpn ) ); Vector vrg( CrossProd( vup, vpn ) ); Ray ray( world->camera->position, vpn + vrg*(x-160)/160 + vup*(y-120)/120 ); But with this method, I get a few problems. 1. I can't conrtol the distance to the viewingplane. 2. How do I find the screen x- and y- coordinates for an arbitary 3d-coordinate. I need this for post-rendering of lens flares. If I always had my camera placed in (0,0,0), looking at (x-160, y-120, 250), there would be no problem. One solution is to, instead of transform the ray with the above algorithm, rotate the world coordinates with vrg,vup,vpn, but this doesn't work when an object is defined by a center and a rotation angle.










-------------------------------------------------------------------------------------------------------- COMPUTER GRAPHICS FAQ Table of Contents --------------------------------------------------------------------------------------------------------- 0. General Information 0.01: Charter of comp.graphics.algorithms 0.02: Are the postings to comp.graphics.algorithms archived? 0.03: How can I get this FAQ? 0.04: What are some must-have books on graphics algorithms? 0.05: Are there any online references? 0.06: Are there other graphics related FAQs? 0.07: Where is all the source?
I have left out Section 0 since it contained no useful algorithms
1. 2D Computations: Points, Segments, Circles, Etc. 1.01: How do I rotate a 2D point? 1.02: How do I find the distance from a point to a line? 1.03: How do I find intersections of 2 2D line segments? 1.04: How do I generate a circle through three points? 1.05: How can the smallest circle enclosing a set of points be found? 1.06: Where can I find graph layout algorithms? 2. 2D Polygon Computations 2.01: How do I find the area of a polygon? 2.02: How can the centroid of a polygon be computed? 2.03: How do I find if a point lies within a polygon? 2.04: How do I find the intersection of two convex polygons? 2.05: How do I do a hidden surface test (backface culling) with 2d points? 2.06: How do I find a single point inside a simple polygon? 2.07: How do I find the orientation of a simple polygon? 3. 2D Image/Pixel Computations 3.01: How do I rotate a bitmap? 3.02: How do I display a 24 bit image in 8 bits? 3.03: How do I fill the area of an arbitrary shape? 3.04: How do I find the 'edges' in a bitmap? 3.05: How do I enlarge/sharpen/fuzz a bitmap? 3.06: How do I map a texture on to a shape? 3.07: How do I detect a 'corner' in a collection of points? 3.08: Where do I get source to display (raster font format)? 3.09: What is morphing/how is it done? 3.10: How do I quickly draw a filled triangle? 3.11: D Noise functions and turbulence in Solid texturing. 3.12: How do I generate realistic sythetic textures? 3.13: How do I convert between color models (RGB, HLS, CMYK, CIE etc)? 4. Curve Computations 4.01: How do I generate a bezier curve that is parallel to another bezier? 4.02: How do I split a bezier at a specific value for t? 4.03: How do I find a t value at a specific point on a bezier? 4.04: How do I fit a bezier curve to a circle? 5. 3D computations 5.01: How do I rotate a 3D point? 5.02: What is ARCBALL and where is the source? 5.03: How do I clip a polygon against a rectangle? 5.04: How do I clip a polygon against another polygon? 5.05: How do I find the intersection of a line and a plane? 5.06: How do I determine the intersection between a ray and a polygon? 5.07: How do I determine the intersection between a ray and a sphere? 5.08: How do I find the intersection of a ray and a bezier surface? 5.09: How do I ray trace caustics? 5.10: What is the marching cubes algorithm? 5.11: What is the status of the patent on the "marching cubes" algorithm? 5.12: How do I do a hidden surface test (backface culling) with 3d points? 5.13: Where can I find algorithms for 3D collision detection? 5.14: How do I perform basic viewing in 3d? 5.15: How do I optimize a 3D polygon mesh? 5.16: How can I perform volume rendering? 5.17: Where can I get the spline description of the famous teapot etc.? 5.18: How can the distance between two lines in space be computed? 5.19: How can I compute the volume of a polyhedron? 6. Geometric Structures and Mathematics 6.01: Where can I get source for Voronoi/Delaunay triangulation? 6.02: Where do I get source for convex hull? 6.03: Where do I get source for halfspace intersection? 6.04: What are barycentric coordinates? 6.05: How do I generate a random point inside a triangle? 6.06: How do I evenly distribute N points on (tesselate) a sphere? 6.07: What are coordinates for the vertices of an icosohedron? 6.08: How do I generate random points on the surface of a sphere? 7. Contributors 7.01: How can you contribute to this FAQ? 7.02: Contributors. Who made this all possible.
I have left out Section 7 also since it contained no useful algorithms
Search e.g. for "Section 6" to find that section. Search e.g. for "Subject 6.04" to find that item. ---------------------------------------------------------------------- Section 1. 2D Computations: Points, Segments, Circles, Etc. ---------------------------------------------------------------------- Subject 1.01: How do I rotate a 2D point? In 2-D, the 2x2 matrix is very simple. If you want to rotate a column vector v by t degrees using matrix M, use M = {{cos t, -sin t}, {sin t, cos t}} in M*v. If you have a row vector, use the transpose of M (turn rows into columns and vice versa). If you want to combine rotations, in 2-D you can just add their angles, but in higher dimensions you must multiply their matrices. ---------------------------------------------------------------------- Subject 1.02: How do I find the distance from a point to a line? Let the point be C (XC,YC) and the line be AB (XA,YA) to (XB,YB). The length of the line segment AB is L: L=((XB-XA)**2+(YB-YA)**2)**0.5 and (YA-YC)(YA-YB)-(XA-XC)(XB-XA) r = ----------------------------- L**2 (YA-YC)(XB-XA)-(XA-XC)(YB-YA) s = ----------------------------- L**2 Let I be the point of perpendicular projection of C onto AB, the XI=XA+r(XB-XA) YI=YA+r(YB-YA) Distance from A to I = r*L Distance from C to I = s*L If r<0 I is on backward extension of AB If r>1 I is on ahead extension of AB If 0<=r<=1 I is on AB If s<0 C is left of AB (you can just check the numerator) If s>0 C is right of AB If s=0 C is on AB ---------------------------------------------------------------------- Subject 1.03: How do I find intersections of 2 2D line segments? This problem can be extremely easy or extremely difficult depends on your applications. If all you want is the intersection point, the following should work: Let A,B,C,D be 2-space position vectors. Then the directed line segments AB & CD are given by: AB=A+r(B-A), r in [0,1] CD=C+s(D-C), s in [0,1] If AB & CD intersect, then A+r(B-A)=C+s(D-C), or XA+r(XB-XA)=XC+s(XD-XC) YA+r(YB-YA)=YC+s(YD-YC) for some r,s in [0,1] Solving the above for r and s yields (YA-YC)(XD-XC)-(XA-XC)(YD-YC) r = ----------------------------- (eqn 1) (XB-XA)(YD-YC)-(YB-YA)(XD-XC) (YA-YC)(XB-XA)-(XA-XC)(YB-YA) s = ----------------------------- (eqn 2) (XB-XA)(YD-YC)-(YB-YA)(XD-XC) Let I be the position vector of the intersection point, then I=A+r(B-A) or XI=XA+r(XB-XA) YI=YA+r(YB-YA) By examining the values of r & s, you can also determine some other limiting conditions: If 0<=r<=1 & 0<=s<=1, intersection exists r<0 or r>1 or s<0 or s>1 line segments do not intersect If the denominator in eqn 1 is zero, AB & CD are parallel If the numerator in eqn 1 is also zero, AB & CD are coincident If the intersection point of the 2 lines are needed (lines in this context mean infinite lines) regardless whether the two line segments intersect, then If r>1, I is located on extension of AB If r<0, I is located on extension of BA If s>1, I is located on extension of CD If s<0, I is located on extension of DC Also note that the denominators of eqn 1 & 2 are identical. References: [O'Rourke] pp. 249-51 [Gems III] pp. 199-202 "Faster Line Segment Intersection," ---------------------------------------------------------------------- Subject 1.04: How do I generate a circle through three points? Let the three given points be a, b, c. Use _0 and _1 to represent x and y coordinates. The coordinates of the center p=(p_0,p_1) of the circle determined by a, b, and c are: A = b_0 - a_0; B = b_1 - a_1; C = c_0 - a_0; D = c_1 - a_1; E = A*(a_0 + b_0) + B*(a_1 + b_1); F = C*(a_0 + c_0) + D*(a_1 + c_1); G = 2.0*(A*(c_1 - b_1)-B*(c_0 - b_0)); p_0 = (D*E - B*F) / G; p_1 = (A*F - C*E) / G; If G is zero then the three points are collinear and no finite-radius circle through them exists. Otherwise, the radius of the circle is: r^2 = (a_0 - p_0)^2 + (a_1 - p_1)^2 Reference: [O' Rourke] p. 201. Simplified by Jim Ward. ---------------------------------------------------------------------- Subject 1.05: How can the smallest circle enclosing a set of points be found? This circle is often called the minimum spanning circle. It can be computed in O(n log n) time for n points. The center lies on the furthest-point Voronoi diagram. Computing the diagram constrains the search for the center. Constructing the diagram can be accomplished by a 3D convex hull algorithm; for that connection, see, e.g., [O'Rourke, p.195ff]. For a direct algorithm, see: S. Skyum, "A simple algorithm for computing the smallest enclosing circle" Inform. Process. Lett. 37 (1991) 121--125. ---------------------------------------------------------------------- Subject 1.06: Where can I find graph layout algorithms? See the paper by Eades and Tamassia, "Algorithms for Drawing Graphs: An Annotated Bibliography," Tech Rep CS-89-09, Dept. of CS, Brown Univ, Feb. 1989. A revised version of the annotated bibliography on graph drawing algorithms by Giuseppe Di Battista, Peter Eades, and Roberto Tamassia is now available from ftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.tex.gz and ftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.ps.gz ---------------------------------------------------------------------- Section 2. 2D Polygon Computations ---------------------------------------------------------------------- Subject 2.01: How do I find the area of a polygon? The signed area can be computed in linear time by a simple sum. The key formula is this: If the coordinates of vertex v_i are x_i and y_i, twice the signed area of a polygon is given by 2 A( P ) = sum_{i=0}^{n-1} (x_i y_{i+1} - y_i x_{i+1}). Here n is the number of vertices of the polygon. References: [O' Rourke] pp. 18-27; [Gems II] pp. 5-6: "The Area of a Simple Polygon", Jon Rokne. To find the area of a planar polygon not in the x-y plane, use: 2 A(P) = abs(N . (sum_{i=0}^{n-1} (v_i x v_{i+1}))) where N is a unit vector normal to the plane. The `.' represents the dot product operator, the `x' represents the cross product operator, and abs() is the absolute value function. [Gems II] pp. 170-171: "Area of Planar Polygons and Volume of Polyhedra", Ronald N. Goldman. ---------------------------------------------------------------------- Subject 2.02: How can the centroid of a polygon be computed? The centroid (a.k.a. the center of mass, or center of gravity) of a polygon can be computed as the weighted sum of the centroids of a partition of the polygon into triangles. The centroid of a triangle is simply the average of its three vertices, i.e., it has coordinates (x1 + x2 + x3)/3 and (y1 + y2 + y3)/3. This suggests first triangulating the polygon, then forming a sum of the centroids of each triangle, weighted by the area of each triangle, the whole sum normalized by the total polygon area. This indeed works, but there is a simpler method: the triangulation need not be a partition, but rather can use positively and negatively oriented triangles (with positive and negative areas), as is used when computing the area of a polygon. This leads to a very simple algorithm for computing the centroid, based on a sum of triangle centroids weighted with their signed area. The triangles can be taken to be those formed by one fixed vertex v0 of the polygon, and the two endpoints of consecutive edges of the polygon: (v1,v2), (v2,v3), etc. The area of a triangle with vertices a, b, c is half of this expression: (b[X] - a[X]) * (c[Y] - a[Y]) - (c[X] - a[X]) * (b[Y] - a[Y]); Code available at ftp://grendel.csc.smith.edu/pub/code/centroid.c (3K). Reference: [Gems IV] pp.3-6; also includes code. ---------------------------------------------------------------------- Subject 2.03: How do I find if a point lies within a polygon? The definitive reference is "Point in Polyon Strategies" by Eric Haines [Gems IV] pp. 24-46. The code in the Sedgewick book Algorithms (2nd Edition, p.354) is incorrect. The essence of the ray-crossing method is as follows. Think of standing inside a field with a fence representing the polygon. Then walk north. If you have to jump the fence you know you are now outside the poly. If you have to cross again you know you are now inside again; i.e., if you were inside the field to start with, the total number of fence jumps you would make will be odd, whereas if you were ouside the jumps will be even. The code below is from Wm. Randolph Franklin with some minor modifications for speed. It returns 1 for strictly interior points, 0 for strictly exterior, and 0 or 1 for points on the boundary. The boundary behavior is complex but determined; in particular, for a partition of a region into polygons, each point is "in" exactly one polygon. See the references below for more detail. int pnpoly(int npol, float *xp, float *yp, float x, float y) { int i, j, c = 0; for (i = 0, j = npol-1; i < npol; j = i++) { if ((((yp[i]<=y) && (y